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Abstract. We describe an algorithm to optimally extract individual spectra of blended sources from a long slit 
spectrum. A semi-analytic model for the spatial profile is used: a Voigt profile for the undersampled core with 
a numerical correction applied to the wings. This is derived from an isolated source which must also be placed 
on the slit, x 2 fitting is used to separate the blended components with maximum weight given to the best data. 
We demonstrate the successful application of the algorithm to data on two X-ray binaries in crowded fields: 
XTE J2012+381 and V404 Cyg. 
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1. Introduction 

A common situation in observational astronomy is to 
seek accurate photometric or spectroscopic information 
about sources within crowded fields, i.e. where the near- 
est neighbours are blended with the sources of interest. 
This is a particularly common situation for X-ray binaries, 
most of which are faint Galactic plane objects. A solution 
to the photometric problem has long been available in 
various i mplem entations of Stetson's DAOPHOT program 
(Stetson 1987 ). For crowded field spectroscopy, however, 
no widely available solution exists. A possible approach to 
the problem is described here and has much in common 
with daophot. This spectroscopic solution also seeks to 
be optimal, in the sense that it should yield the highest 
signal-to-noise spectra available from the data. It thus de- 
rives from optimal extraction methods developed for un- 
blended spectra (Home 1986 ; Marsh 198£ ) . We developed 
the method initially in the context of an observation of the 



X-ray binary XTE J2012+381 (Hynes et al. |1999j. It has 
also provided useful for another X-ray binary, V404 Cyg 



(Hynes et al. 2002). We illustrate the application of the 
method to both these sources in Section [| 

Since developing the method, another approach has 



been presented by Buie & Grundy (2000). This is also 
an optimal extraction, but differs in that no template 
source is used; instead the spatial profile is reconstructed 
iteratively from the blended sources. These authors also 
use a purely numerical profile, rather than the semi- 
analytic model we adopt. Their approach was developed 
for ffiST/NICMOS observations of a Pluto and Charon; 
hence these requirements were necessary. For ground- 
based observations of point sources in crowded fields, us- 



ing a template source and a semi-analytic model profile 
should allow our method to be applied to poorer quality 
data. We should emphasise that the assumptions made 
here do only apply to point sources; this method cannot 
be used for sources with extended or resolved structure. 

2. Preliminary Image Processing 

Before applying the method, the images should have been 
debiased and flat-fielded using standard techniques. It is 
assumed hereafter that the image is oriented with the dis- 
persion direction roughly horizontal and the slit is aligned 
with columns. 

Sky subtraction should also have been performed and 
to ensure correct weighting across the profile, this should 
be done in a way which yields a image of the fitted sky 



background. Home (1986) discusses these issues and we 
essentially follow his methods. Our implementation fits a 
low-order polynomial to each column with Home's itera- 
tive cosmic ray rejection. Regions containing sources are 
masked out. We then retain both the subtracted image 
and an image containing the fitted sky. 

It is also necessary to do wavelength calibration sep- 
arately using normal procedures. The method employed 
for the spectra presented in Section ^ was to extract and 
wavelength calibrate a spectrum of the template source 
using normal IRAF methods, and then assume the same 
wavelength-pixel correspondence for all the objects. This 
is acceptable for these data, since the spectra considered 
do have the image of the slit precisely aligned with the 
CCD columns, but in general a more rigorous procedure 
may be needed. 
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3. The Spatial Profile 



3.2. Fitting the profile 



3.1. The profile model 

The first step in deblending spatial profiles is to define how 
the profile of a single point source should look. Several 
metho ds hav e been proposed for this for a single source. 
Home ( 1986 ) normalises the observed count rates for each 
column, then fits a low order polynomial to each row of 
the data. This has the advantage of making no a priori 
assumptions about the spatial profile, but is inappropri- 
ate for spectrum deblending because his spatial profile is 
only defined in pixel steps; it cannot easily be transferred 
to another source which samples the spatial profile differ- 
ently. In principle, one could interpolate between pixels, 
but this is dangerous when the spatial profile is under- 
sampled, as is often the case (even where the seeing is 
poor enough that the spatial profile would become over- 
sampled, it is common to bin to reduce readout noise). 
Home's method is also only appropriate for spectra with 
small distortion, i.e. nearly parallel to CCD rows. Marsh 
( |1989D describes an extension of this empirical approach 
which works even for very distorted spectra. His method 
also produces well sampled spatial profiles which can be 
transferred to another source on the slit, taking advan- 
tage of the fact that strong distortion of the spectrum will 
lead to different columns sampling the spatial profile dif- 
ferently. This could be used as the basis for a deblending 
algorithm. In some cases, however, the distortion is very 
slight so this method will not work either; even by combin- 
ing columns, the spatial profile is not well sampled. This 
was the case for the two datasets discussed in Section [|, 
and so led us to develop a different method. 

The alternative is simpler to visualise, since it involves 
fitting a analytical profile in the spatial direction rather 
than a polynomial in the dispersion direct ion. Analytic 
spatial profiles have been applied by Urry ( 1988 ) to IUE 
data using a Gaussian profile. For ground based data a 
Gaussian model will often give a good fit to the seeing 
dominated cores of spatial profiles, but will usually un- 
derpredict the extended wings present due to instrumen- 
tal imperfections. An effective refinement to this is to use 
a Voigt profile, a convolution of a Gaussian profile with a 
dispersion profile (e.g. Gray 1992): 

exp( — uf) 



Pi(ui, a) oc 



r dui. 



(1) 



where Pj is the profile value, Uj = Xi/w, Xi is the pixel po- 
sition, w is the profile width parameter and a is the Voigt 
damping parameter. This function provides the desired 
Gaussian cores, but also gives an independently variable 
extension in the wings, controlled by the Voigt damping 
parameter which ranges from (a pure Gaussian profile) 
to 1 (a Lorentzian profile). It was found that a Voigt profile 
provided a very good fit to the spatial profiles considered 
in Section ^. The resulting systematic residuals to the fit 
are very small; a numerical correction method similar to 
that used by DAOPHOT can be used to refine the fit even 
more as described in the following subsection. 



The first step in the fitting process is to perform an un- 
constrained fit to each column of the template data, with 
the centre, Gaussian width, Voigt damping parameter and 
profile scaling as free parameters. We have used a downhill 
simplex (amoeba) algorithm to obtain the best-fit param- 
eters (e.g. Press et al. 1992). No attempt is made to reject 
bad pixels or cosmic rays at this stage, as these will show 
up as anomalous parameter values. Of these parameter 
values, it is the profile centre, width and damping that 
define a normalised profile. All of these parameters can be 
expected to vary smoothly in wavelength, so after obtain- 
ing fitted values, a low order polynomial in wavelength is 
fitted to each parameter, rejecting anomalous values that 
may have been affected by cosmic rays. If the wavelength 
range is small then a zero-order fit (i.e. a constant) may 
suffice. If the data quality is not sufficient to constrain the 
fits to a single row of the image then several rows can be 
binned, since the wavelength dependence is likely to be 
slow, and the derived polynomial resampled back to the 
original resolution. These polynomial fits, together with 
the normalisation such that JV P{ = 1 define the spatial 
profile as a function of wavelength. 

If the difference in brightness of the two sources is very 
large then it becomes crucial that the wing of the brighter 
source be very well fitted to avoid contamination of the 
fainter source. This can be achieved using a numerical cor- 
rection to the (Voigt) model profile. The error involved in 
assuming a Voigt profile is usually only important in the 
extreme wings of the profile: the seeing-dominated core 
is typically well fitted by this model. We can therefore 
define a numerical correction factor (observed profile di- 
vided by model) at the sampled points and interpolate 
this to obtain a general correction function for any pixel 
sampling. Since this involves the extreme wings of the 
line, where signal-to-noise will be poor, such a correction 
typically cannot be defined meaningfully as a function of 
wavelength; instead only a single, averaged, wavelength 
independent correction can be determined. A symmetric 
profile is also assumed. These assumptions appear ade- 
quate for the spectra we have considered; hence the cor- 
rection is simply a function of distance from the profile 
centre. Given higher quality data, an asymmetric, wave- 
length dependent function could readily be used, and this 
generalisation is straightforward. 

Assuming that there are no changes in focus along the 
slit (or at least the restricted part of it which is used) 
and that the image scale (and hence the separation be- 
tween stellar images) is independent of wavelength, the 
only additional information required to define the profiles 
of the blended sources is the separation of each from the 
template source. For spectra with minimal distortion, this 
is easy to obtain by performing an initial fit to the av- 
erage of several rows of data, leaving the centres of the 
two blended profiles as free parameters. Distorted spectra 
would require that the distortion (known from the fits to 
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the template profile) be removed by resampling the data 
for this initial fit only. 



4. Spectral Extraction 

Assuming that the fit to the spatial profile for the tem- 
plate source is satisfactory, the centre position and nor- 
malised profile as a function of wavelength are known, 
and the spectral extraction itself can be performed. This 
is simply a case of finding the scaling of each component, 
although complicated by the need to reject cosmic ray 
events. Home (1986) describes his method in terms of an 
optimally weighted sum of counts over an aperture con- 
taining all of the stellar light. He notes the equivalence 
between this weighted sum method and fitting a known 
profile to data of known variance, and it can be shown that 
his optimally weighted sum yields the solution of minimum 



X ■ Using this equivalence we can generalise the method 
to blended spectra, by minimising the x 2 °f the fit with 
respect to the two profile scalings. This is straightforward 
and for the simple problem of scaling two profiles, the so- 
lution of minimum x can be calculated analytically; see 
Appendix [a|. This is done iteratively, to allow rejection 
of cosmic rays. Given that there are two components to 
be fitted, however, the combined profile is unknown, and 
so more care is needed in cosmic ray rejection than for a 
single source. 



5. Cosmic Ray Rejection 

One of the great strengths of optimal extraction algo- 
rithms is that the spatial profile is known. Cosmic rays 
or bad pixels, even when on top of the spectrum, will dis- 
tort the spatial profile, and hence can be recognised and 
rejected. This means that instead of removing cosmic rays 
from the extracted one-dimensional spectra, and rejecting 
a whole pixel step in wavelength, they can be removed 
from the two-dimensional spectra, retaining the informa- 
tion in the uncontaminated pixels. While this process is 
easy to do by eye, it proves difficult to train a computer 
to recognise cosmic rays automatically, especially when 
multiple profiles are being fit. Given a first approximation 
to the profile, bad values can be rejected by an iterative 
sigma-clipping algorithm (Home 1986 ) . The strategies we 
have used to obtain the crucial first approximation are 
described here. 

Beginning with the simplest case of recognising con- 
tamination of a single profile, we rely on the fact that the 
scaling of the profile can be estimated from a single pixel 
value; hence each pixel gives an independent estimate of 
the integrated flux. A contaminated pixel will then give 
an anomalous estimate. If this estimate is severely wrong, 
as is often the case with cosmic rays, then a straight av- 
erage may be very much in error, leading to subsequent 
rejection of good data. A better method is to take the me- 
dian of the single pixel estimates as the first estimate for 
the integrated flux. This then gives a first approximation 



to the fitted profile which forms the basis for subsequent 
refinement with iterative sigma-clipping. 

Where two independently fitted profiles are involved 
more care is needed. The only case considered here is that 
one source is brighter than the other at all wavelengths 
in the spectrum (true for all data the method has so far 
been applied to.) For the brightest source, a profile region 
is selected from the midpoint between the two sources to 
an arbitrary cutoff on the other side. Within this region, 
the profile of the brighter source is expected to dominate, 
so the median ratio method described above for single 
sources can be applied. This gives an approximate flux 
scaling for the brighter stellar profile, which can then be 
subtracted. The median ratio method is then applied to 
the remaining profile of the fainter source to estimate a 
flux scaling for this. With an approximate scaling for both 
profiles, iterative sigma-clipping can be used to refine the 
list of rejected pixels. 

Clearly if three or more blended sources are involved 
then the situation is more complex, and there are more 
possible permutations of relative brightness to consider. 
This makes it harder to treat generally and specific cases 
will need to be considered individually. 



6. Examples 

6.1. XTE J2012+381 

XTE J2012+381 presents an especially difficult challenge 
for spectroscopy. The X-ray source was originally identi- 
fied with what appeared to be a normal F star. It was only 
upon closer examination of images of the field taken in 
good seeing that it became clear a second fainter star was 
present (Hynes et al. 1999). This was heavily blended with 
the F star, but lay closer to the precise radio counterpart 
and now appears the most likely optical counterpart to 
the X-ray source. We obtained spectroscopy of these two 
stars with the William Herschel Telescope (WHT) on 1998 
July 20, with the slit aligned to pass through both, as well 
as through a third unblended star. These data provided 
the original motivation for development of this method. A 
f uller description of the data is provided by Hynes et al. 
(119991) . We focus here on aspects relevant to the deblcnd- 
ing algorithm. 

The severity of the blending can be clearly seen in Fig. 
||. The suggested counterpart to the X-ray source is the 
fainter of the two blended stars. Clearly to successfully 
extract its spectrum requires a good fit to the wing of the 
brighter star. Fortunately, the signal-to-noise ratio of the 
two-dimensional image was quite good, so it was possible 
to apply the full model described above, i.e. a Voigt pro- 
file with a numerical correction to the wings. The latter 
significantly improved the fit and with this no significant 
discrepancy was seen between the model profile and the 
data. The fits to the Voigt profile parameters are shown 
in Fig. ||. Both the width and damping parameter are rel- 
atively well constrained as a function of wavelength. 
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Fig. 1. Fit to spatial profiles in the red part of the WHT 
spectrum of XTE J2012+381. Data are shown by points, 
the model profile is dashed and the dotted lines indicate 
the De-blended components. Both data and fits have been 
summed over 20 pixels in the dispersion direction to re- 
duce noise and illustrate the quality of fit achieved. On 
the left is the fit to the spatial profile of the template star. 
On the right is the two-profile fit to the blended stars. 



The extracted spectra are shown in Fig. ^. There are 
a shortage of spectral features in both spectra, but the 
extraction does clearly distinguish the much redder spec- 
trum of the proposed X-ray source and the contrast be- 
tween Ha absorption in the F star and possible weak emis- 
sion in the X-ray source is dramatic. The r eality of this 
emission is discussed by Hynes et al. ( 1999| ); suffice here 
to note that no obvious artifacts are present, but that 
these data are perhaps not the best test of the method. 
Fortunately our second target, V404 Cyg, does provide a 
better test, albeit in a less heavily blended situation. 



6.2. V404 Cyg 

Another X-ray binary for which this deblending method is 
appropriate is V404 Cyg, as this has a fainter companion 
star 1.5arcsec away. There is no physical association be- 
tween the two. Observations were obtained of V404 Cyg 
in quiescence using the WHT on 1999 Jul y 6-8. These are 
described in more detail in Hynes et al. ( 2002 ). To sum- 
marise, the goal was to perform spectroscopy near Ha with 
as high a time resolution as possible ( ~ 240 s) to examine 
spectral variability; hence we used a low spectral resolu- 
tion to minimise readout noise and a wide slit to minimise 
slit losses. Figure [| shows two examples of averaged (in 
wavelength) spatial cuts with respect to the centre of the 
template star. The images with the best and worst seeing 
from the night of July 7/8 were chosen. The fainter star 
is clearly blended with V404 Cyg even in the best image 
(seeing ~ 0.8arcsec.) 

The optimal deblending algorithm was applied with 
a simplification. Because the goal of this project was to 
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Fig. 2. Wavelength dependence of the Voigt profile fit- 
ting parameters, the Gaussian width and the Voigt damp- 
ing parameter. Points mark fits to individual pixels in the 
dispersion direction, with no binning. The solid lines are 
low-order polynomial fits to the parameters. 

obtain the highest time-resolution practical, the signal- 
to-noise ratio of individual frames was not high. Because 
of this, it was not possible to perform a wavelength de- 
pendent fit and so a single average profile (and a profile 
correction) were calculated. Given the limited wavelength 
coverage of these data this is an acceptable approximation. 

Figure || shows the average extracted spectra of the two 
blended stars from July 7/8. The spectra appear to have 
been separated cleanly and there is no obvious crosstalk 
around Ha between the double-peaked disc emission in 
V404 Cyg and the narrower (stellar) absorption line in the 
blended star. There is also a lot of structure in the spec- 
trum of V404 Cyg (e.g. 6600-6900 A) which is not present 
in the spectrum of the blended star which appears largely 
featureless. As V404 Cyg is brighter, the signal-to-noise 
ratio of its spectrum should be higher and hence these fea- 
tures must be real; they likely arise from the KOIV com- 
panion star which dominates in quiescence. The blended 
star, in contrast, appears to be of earlier spectral type, 
probably F (Hynes et al. |2002| )B 

Figure ^| shows continuum light curves from both stars 
for the night of July 7/8. Spectra were scaled relative to 
the template star to approximately remove variable slit 
loss and transparency effects. Clearly the light curve of 
the blended star shows little or no contamination from 



1 To the author's knowledge, no spectrum of this star has 
previously been isolated. 



R. I. Hynes: An optimal extraction of spatially blended spectra 



5 



7 r 




Ej , , I , , , , I , , , , I , , , , I , , , , I , , lZ 

6000 6500 7000 7500 8000 

Wavelength (A) 

Fig. 3. WHT spectra of the F star and the faint 
red companion believed to be the optical counterpart 
of XTE J20 12+381. Both have been binned x4 in 
wavelength for clarity. Atmospheric absorption features 
(marked ©) have been corrected for in both spectra. The 
only prominent line in either spectrum is Ha: absorption 
in the brighter star and possible weak emission in the sug- 
gested counterpart to the X-ray source. 




L , , I , , , I , , , I , , , I , , , I , , , I , , , L_ 

6000 6200 6400 6600 6800 7000 7200 7400 

A (A) 

Fig. 5. Deblended spectra of V404 Cyg and the blended 
star. All the spectra from 1999 July 7/8 have been av- 
eraged. The deblcnding algorithm has appears to have 
separated the spectra cleanly with the Ha emission from 
V404 Cyg clearly distinguished from absorption in the 
blended star. The absorption feature in the blended star 
is not exactly aligned with the dip in the line profile of 
V404 Cyg. 




Pixel offset from template star 

Fig. 4. Best (solid) and worst (dashed) average spatial 
profiles of V404 Cyg from 1999 July 7/8. These are the 
observed profiles not models. Columns containing cosmic 
rays have been excluded from the average. Peaks are from 
left to right, the blended companion star, V404 Cyg and 
the profile template star. The profiles have been scaled 
such that the peak flux from the template star is unity. 

the large flare occurring at the beginning of the night in 
V404 Cyg, indicating again that there is no detectable 
cross-talk between the two spectra. 

7. Conclusions 

We have demonstrated an algorithm for extracting sepa- 
rate spectra of sources whose spatial profiles are blended 
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Fig. 6. Continuum light curves of V404 Cyg and its 
blended companion. 

in a long slit spectrum given an isolated template source 
with which to define the spatial profile. We implement 
this method using x 2 fitting with rejection of cosmic rays, 
resulting in a generalisation of the optimal extraction 
method to the multi-object case. Finally we have demon- 
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strated the application of the method to two X-ray bina- 
ries in crowded fields that require separation from nearby 
contaminating sources. 

The method as presented here has some limitations: 

1. A very simple method of wavelength calibration is 
used, i.e. assuming that the slit is exactly aligned 
with CCD columns. A more rigorous treatment would 
be able to use a two-dimensional wavelength solution 
when necessary. 

2. The assumption of a Gaussian core to the profile may 
not be valid if the image size is determined by poor 
focus rather than seeing, or if the exposure time is 
short enough 1 s) that the seeing does not average 
to a Gaussian form. 

3. At present the method is targeted at two blended pro- 
files. A generalisation to three or more sources would 
not only require an extension of Appendix [A], but de- 
velopment of a method for identifying cosmic rays su- 
perposed on such a complex profile. 

An implementation of the method described using IDL 
routines operating on FITS images is available from the 
author. 
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reduces to finding the flux scaling of each profile. The 
optimal solution is equivalent to minimising \ 2 , with care 
to reject any cosmic rays on the profile. 

Consider N model profiles (J = 1 . . . N) to be fitted 
to an observed profile Di (variance Vi), yielding N scalings 
Fj. The badness of fit statistic is then: 



E 



(A - J2j F j p jj) 

Vi 



(AT) 



and the optimal solution requires dx 2 /dFj = for all j. 
In general, these N constraints will yield N simultaneous 
equations with N unknowns: F\ . . . Fj\f. This can be solved 
as a straightforward problem in linear algebra without re- 
quiring a numerical minimisation. Consider the two-profile 
case. For convenience write F\ = /, F% = g, P^i = Pi and 
Pi 2 = qi- Define the following statistics: 

A = E^ B = E tiK c = X>*/ y < 



D 



Y.PiDi/Vi E = ^q i D i /V i . 



(A.2) 



The two simultaneous equations then take the form 

Af -D + Cg = Bg-E + Cf = (A.3) 

with solution 



, D/C-E/B E/C-D/A 
f = -CTF, 7777T 9~ 



(A.4) 



A/C-C/B * BIG -Cj A 
Formal errors can also be straightforwardly propagated 



from Eqn. A.3 or A.4 



^ \ dDi 



E 



dg 
8D, 



V,. = 
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AB-C 2 
A 

AB-C 2 



(A.5) 



Appendix A: Analytic \ 2 fitting of multiple 
profiles 

For the simple case considered here, the spatial profile of 
blended sources is known a priori from a template source 
and it is assumed that the offsets of these blended sources 
relative to the template are known. The extraction then 



